SSA

table <- read.csv2("ssa.csv")
long <- table$COL1
short <- na.omit(table$COL7)

Long series

plot(long, type='l')

spec.pgram(long, detrend = TRUE, fast = FALSE, log='no', taper=0)

Предполагем, что достаточно сгруппировать 8 компонент, т.к. возможно, ранг равен 8 (линейный тренд + 3 косинуса).

Как выбирать L

Так как точной разделимости в реальных рядах не бывает, нас интересует асимтотическая разделимость.

В большинстве случаев скорость асимптотической разделимиости пропорциональна \(1/\min(L, K)\), \(K = N - L + 1\), поэтому выбираем \(L\) близкое к \(N/2\).

Если мы хотим отделять косинус с некоторым периодом, то для улучшения разделимости рекомендуется брать окно кратное периоду (следует из условий точной слабой разделимости).

Слабая и сильная разделимость

• Слабая разделимость – ортогональность траекторных подпространств
• Сильная азделимость – слабая разделимость + в разных компонентах (SVD) разные собственные числа

Decomposition

ssa.long <- ssa(long, kind='1d-ssa')

Собственные числа:

\(\lambda_i = \|\mathbf{X}_i\|^2_{\mathcal{F}}\), \(\lambda_i/ \|\mathbf{X}\|^2_{\mathcal{F}}\) – вклад компоненты

plot(ssa.long)

Как идентифицировать тренд?

Хотим выделить все элементарные компоненты, которые имеют медленно меняющиеся собственные вектора.

Собственные ветора образуют базис траекторного подпространства, следовательно, они имеют вид соответствующих компонент ряда.

Собственные вектора:

plot(ssa.long, type='vectors', idx = 1:15)

Первые два – линейный тренд, дальше идут периодические компоненты.

Как идентифицировать периодичность

Посмотрим на скаттерплоты пар собственных векторов.

• Знаем, что экспоненциально модулированный косинус соответствует двум собственным тройкам с примерно одинаковыми собственными числами (\(Lw\) – целые). Отсюда следует, что эти SVD компоненты скорее всего идут друг за другом
• Эти две периодические компоненты (собственные вектора), отличаются на \(\pi/2\) по фазе (при целых \(Lw\)).

Таким образом, видимые Т-угольники соответствуют периодическим компонентам с периодом Т.

plot(ssa.long, type='paired', idx= c(1:20))

Заметны компоненты с периодами 20, 10 и 6. Есть компонента, похожая на период 2.4, компонента похожая на период ~6.

Элементарные восстановленные компоненты

plot(ssa.long, type='series', groups = list(1:2, 3:4, 5:6, 7:8, 9:10, 11:12))

Оценка параметров (ESPRIT)

\(x_n = \sum\limits_{j = 1}^r C_j \mu_j^n + \epsilon_n\) – ряд, \(\mu_j = \rho_je^{2\pi i \omega_j}\) – корни х.п.

Алгоритм ESPRIT-LS:

• Применив SSA и получив SVD разложение, составляем r собственных векторов в матрицу \(\mathbb{U}\)
• \(\widehat{\mathbb{D}} = \underline{\mathbb{U}}^-\bar{\mathbb{U}}\)
• Собственные чилса – оценки сигнальных корней

parestimate(ssa.long, groups = list(1:12), method = "esprit-ls")

   period     rate   |    Mod     Arg  |     Re        Im
   10.008   0.000360 |  1.00036   0.63 |  0.80960   0.58760
  -10.008   0.000360 |  1.00036  -0.63 |  0.80960  -0.58760
   20.006   0.000299 |  1.00030   0.31 |  0.95137   0.30903
  -20.006   0.000299 |  1.00030  -0.31 |  0.95137  -0.30903
      Inf   0.000104 |  1.00010   0.00 |  1.00010   0.00000
      Inf  -0.000109 |  0.99989   0.00 |  0.99989   0.00000
    5.999  -0.000352 |  0.99965   1.05 |  0.49963   0.86583
   -5.999  -0.000352 |  0.99965  -1.05 |  0.49963  -0.86583
    2.217  -0.000564 |  0.99944   2.83 | -0.95263   0.30227
   -2.217  -0.000564 |  0.99944  -2.83 | -0.95263  -0.30227
    6.665  -0.001292 |  0.99871   0.94 |  0.58685   0.80810
   -6.665  -0.001292 |  0.99871  -0.94 |  0.58685  -0.80810

Компоненты смешались. Как понять, это слабая или сильная разделимость?

Нет слабой разделимости – отсутствие ортогональности траекторных подпространств.

• Посмотреть на матрицу взвешенных корреляций \(\rho_{ij}(X_N^{(i)},X_N^{(j)}) = \frac{(X^{(i)}_N, X^{(j)}_N)_w}{\|X^{(i)}_N\|_w\|X^{(j)}_N\|_w}\).

По такой матрице поймем, хорошо ли мы разделили компоненты. Поймем, какие компоненты делить не стоит.

• Можно посмотреть на частотные коэффициенты корреляции (достаточное условие слабой отделимости)

W-корреляционная матрица

plot(wcor(ssa.long, groups = as.list(1:15)))

Матрица практически диагональная – тренды и сезонные компоненты практически w-ортогональны

Нет сильной разделимости – одинаковые собственные числа оказались в разных компонентах

• Можно посмотреть на вклады компонент, ступеньки будут соответствовать близким собственным числам

Reconstruction

В предположении, что ранг ряда 8, получим

long.reconstructed <- reconstruct(ssa.long, list(
  trend = c(1, 2),
  periodic = 3:8))
res.long <- residuals(long.reconstructed)
par(mfrow=c(4,1))
plot(long, type = 'l')
plot(long.reconstructed$trend, type='l', ylab = "trend")

plot(long.reconstructed$periodic, type='l', ylab = "periodic")
plot(res.long, type='l')

Автоковариационная функция

acf(res.long, lag.max = 50)

Шум похож на белый.

Корни ХП

Min-norm ЛРФ

• \(P_i\) – столбцы траекторной матрицы
• \(\underline{P}_{i}\) – вектор \(P_i\) без последней координаты
• \(\pi_i\) – последняя координата \(P_i\)

Коэффициенты min-norm ЛРФ:

\(\mathcal{R} = \frac{1}{1-\nu^2}\sum\limits_{i=1}^{d}{\pi_i \underline{P}_{i}}\), где \(\nu^2 = \pi_1^2 + … + \pi_d^2\)

\(x_{i + d} = \sum\limits_{k = 1}^da_kx_{i + d - k}\), \(a_d \neq 0\) – минимальная ЛРФ

\(P(\mu) = \mu^d - \sum\limits_{k = 1}^da_k\mu^{d - k}\) – характеристический полином ЛРФ

\(\{\mu_j\}_{j = 1}^p\) – корни х.п. кратности \(k_j\), \(\sum k_j = d\), тогда

\(x_n = \sum\limits_{m = 1}^p(\sum\limits_{j = 0}^{km - 1} C_{m_j} n^j)\mu_m^n\)

min-norm ЛРФ и ее корни (предполагаем, что ранг ряда 8)

lr <- lrr(ssa.long, groups = list(1:8))
plot(lr)

r <- roots(lr)

Хотим построить формулу для сигнала ряда с шумом.

Взяв большое окно и построив min-norm ЛРФ, найдем много корней х.п.

Как выбрать из них те r корней, что определяют структуру сигнала?

Найдем сигнальные корни, воспользовавшись утверждением, что у min-norm ЛРФ все лишние корни по модулю меньше единицы.

mods <- Mod(r)
print(r[mods >= 1])

[1] 0.8098255+0.5875445i 0.8098255-0.5875445i 0.9513031+0.3091901i 0.9513031-0.3091901i 1.0000664+0.0000000i

Перевернем ряд

lr.back <- lrr(ssa.long, groups = list(1:8), reverse = TRUE)
r.back <- roots(lr.back)
mods.back <- Mod(r.back)
print(r.back[mods.back >= 1])

[1] 1.000001+0i

Модули главных корней

main.roots <- c(r[Mod(r) >= 1], 1/r.back[Mod(r.back) >= 1])
print(Mod(main.roots))

[1] 1.0005128 1.0005128 1.0002881 1.0002881 1.0000664 0.9999993

Соответствующие периоды

print(2*pi/Arg(main.roots))

[1]  10.01067 -10.01067  19.99437 -19.99437       Inf       Inf

plot(main.roots, xlim = c(-1,1), ylim = c(-1,1), asp = 1, xlab = "Re", ylab = "Im")
draw.circle(0,0,1,border = 'red')

Нашлось только 6 сигнальных корней, а должно быть 8.

Попробуем другой способ:

Составим систему и найдем коэффициенты по МНК

complex.LS <- function(series, series.roots){
   X <- mapply(FUN = function(root){root^(1:length(series))}, series.roots)
   Y <- cbind(series)
   list(coef = solve(Conj(t(X))%*%X) %*% Conj(t(X)) %*% Y,
        X = X)
}
coeffs.estimates <- complex.LS(long.reconstructed$trend + long.reconstructed$periodic, r)

Воспользуемся тем, что перед сигнальными корнями должны стоять максимальные коэффициенты

idx.roots <-order(Mod(coeffs.estimates$coef), decreasing = TRUE)[1:8]
main.roots<-r[idx.roots]
res.coeffs <- coeffs.estimates$coef[idx.roots]

Модули главных корней

print(Mod(main.roots))

[1] 0.9999324 1.0000664 0.6414650 1.0002881 1.0002881 0.9865705 0.9865705 0.9874613

Соответствующие периоды

print(2*pi/Arg(main.roots))

[1]       Inf       Inf       Inf -19.99437  19.99437  13.87755 -13.87755  13.58822

plot(main.roots, xlim = c(-1,1), ylim = c(-1,1), asp = 1, xlab = "Re", ylab = "Im")
draw.circle(0,0,1,border = 'red')

relation <- function(n){
  rbind(res.coeffs) %*% t(mapply(FUN = function(root){root^n}, main.roots))
}

plot(long.reconstructed$trend + long.reconstructed$periodic, type = 'l', ylab = "series")
lines(x = 1:1000, y = Re(relation(1:1000)), col = 'red')

Short series

plot(short, type='l')

spec.pgram(short, detrend = TRUE, fast = FALSE, log='no', taper=0, xaxt = 'n')
l <- seq(0, 0.5, 0.05)
axis(1, at = l, labels = l)

Ожидаем, что ранг ряда как минимум 4 (тренд + синус).

Decomposition

short.ssa <- ssa(short, kind='1d-ssa')

Собственные числа

plot(short.ssa)

Собственные вектора:

plot(short.ssa, type='vectors', idx = 1:20)

plot(short.ssa, type='paired', idx= c(1:20))

parestimate(short.ssa, groups = list(3:4, 5:6, 7:8, 11:12, 13:14))

$F1
   period     rate   |    Mod     Arg  |     Re        Im
   17.587   0.000000 |  1.00000   0.36 |  0.93686   0.34971

$F2
   period     rate   |    Mod     Arg  |     Re        Im
    5.151   0.000000 |  1.00000   1.22 |  0.34386   0.93902

$F3
   period     rate   |    Mod     Arg  |     Re        Im
    6.265   0.000000 |  1.00000   1.00 |  0.53780   0.84307

$F4
   period     rate   |    Mod     Arg  |     Re        Im
    3.133   0.000000 |  1.00000   2.01 | -0.42103   0.90705

$F5
   period     rate   |    Mod     Arg  |     Re        Im
    2.363  -0.000000 |  1.00000   2.66 | -0.88605   0.46360

plot(wcor(short.ssa, groups = as.list(1:20)))

short.reconstructed <- reconstruct(short.ssa, list(
  trend = c(1, 2),
  periodic = c(3:8,11:14)))
res.short <- residuals(short.reconstructed)
par(mfrow=c(4,1))
plot(short, type = 'l')
plot(short.reconstructed$trend, type='l', ylab = "trend")

plot(short.reconstructed$periodic, type='l', ylab = "periodic")
plot(res.short, type='l')

acf(res.short, lag.max = 50)

Projection SSA

Двойное центрирование

\(\mathcal{A}(\mathbf{X}) = [\mathcal{E}_1(\mathbf{X}) : \ldots : \mathcal{E}_1(\mathbf{X})]\), где \(\mathcal{E}_1(\mathbf{X})\) – вектор из средних по строкам

\(\mathcal{B}(\mathbf{X}) = [\mathcal{E}_{12}(\mathbf{X}) : \ldots : \mathcal{E}_{12}(\mathbf{X})]\), где \(\mathcal{E}_{12}(\mathbf{X})\) – вектор, в котором j-ая компонента равна среднему компонент \(X_j^{(c)} = X_i - \mathcal{E}_1(\mathbf{X})\).

При двойном центрировании применяем SVD к матрице \(\mathbf{A} = \mathcal{A}(\mathbf{X}) + \mathcal{B}(\mathbf{X})\)

Получим разложение \(\mathbf{X} = \mathbf{A} + \sum\limits_{i = 1}^d\sqrt{\lambda_i}U_iV_i^\mathrm{T}\), где \((\lambda_i, U_i, V_i)\) соответствуют \(S^* = (\mathbf{X} - \mathbf{A})(\mathbf{X} - \mathbf{A})^\mathrm{T}\).

Так как исходный тренд линейный, двойное центрирование поможет отделить линейную часть от колебаний (линейная часть попадет в \(\mathbf{A}\)).

short.projection.ssa <- ssa(x = short, column.projector = 'centering', row.projector = 'centering')

plot(short.projection.ssa)

plot(short.projection.ssa, type='vectors', idx = 1:20)

По виду собственных векторов видно, что тренд отделился лучше.

plot(short.projection.ssa, type='paired', idx= c(1:20))

parestimate(short.projection.ssa, groups = list(3:4, 5:6, 7:8, 11:12, 13:14), method ="esprit-ls")

$F1
   period     rate   |    Mod     Arg  |     Re        Im
   19.760  -0.003204 |  0.99680   0.32 |  0.94683   0.31165
  -19.760  -0.003204 |  0.99680  -0.32 |  0.94683  -0.31165

$F2
   period     rate   |    Mod     Arg  |     Re        Im
    5.332  -0.021225 |  0.97900   1.18 |  0.37444   0.90456
   -5.332  -0.021225 |  0.97900  -1.18 |  0.37444  -0.90456

$F3
   period     rate   |    Mod     Arg  |     Re        Im
    6.238  -0.037821 |  0.96289   1.01 |  0.51434   0.81400
   -6.238  -0.037821 |  0.96289  -1.01 |  0.51434  -0.81400

$F4
   period     rate   |    Mod     Arg  |     Re        Im
    3.194  -0.016631 |  0.98351   1.97 | -0.37971   0.90725
   -3.194  -0.016631 |  0.98351  -1.97 | -0.37971  -0.90725

$F5
   period     rate   |    Mod     Arg  |     Re        Im
    2.359  -0.067364 |  0.93485   2.66 | -0.83013   0.42993
   -2.359  -0.067364 |  0.93485  -2.66 | -0.83013  -0.42993

short.reconstructed <- reconstruct(short.projection.ssa, list(
  trend = c(1, 2),
  periodic = c(3:8,11:14)))
res.projection <- residuals(short.reconstructed)
par(mfrow=c(4,1))
plot(short, type = 'l')
plot(short.reconstructed$trend, type='l', ylab = "trend")

plot(short.reconstructed$periodic, type='l', ylab = "periodic")
plot(res.projection, type='l')

acf(res.projection, lag.max = 50)

Стало немного лучше, но все-таки что-то осталось.

Toeplitz SSA vs Basic SSA

Для коротких рядов, которые предполагаются стационарными рекомендуется заменять матрицу \(\mathbf{S} = \mathbf{X}\mathbf{X}^\mathrm{T}\) на автоковариационную матрицу \(\widetilde{\mathbf{C}}\) на этапе разложения

\(\widetilde{c}_{ij} = \frac{1}{N - |i - j|} \sum\limits_{m = 1}^{N - |i - j|} x_m x_{m + |i - j|}\), \(1 \leq i, j \leq L\).

Замечание:

Полученное раложение \(\mathbb{X} = \sum\limits_{i = 1}^L\sigma_iP_iQ_i^\mathrm{T}\), где \(P_i\) – с.в. \(\widetilde{\mathbf{C}}\), образуют ортонормированный базис в \(\mathbb{R}\), \(Q_i = \frac{\mathbb{X}^\mathrm{T}P_i}{\|\mathbb{X}^\mathrm{T}P_i\|}\). Это разложение не является сингулярным разложением, \(\sigma_i = \|\mathbb{X}^\mathrm{T}P_i\|\) вообще говоря не собственные числа \(\widetilde{\mathbf{C}}\).

Пример: ряд \(x_n = e^{\alpha n} cos(2\pi n/7) + \epsilon_n\)

set.seed(100)
eps <- rnorm(100)
alpha <- seq(0.0001, 0.01, 0.0001)
x <- mapply(FUN = function(a) {exp(a * (1:100)) * cos(2*pi*(1:100)/7) + eps}, alpha)
basic <- apply(x, 2, FUN = function(series){
  s <- ssa(series)
  rec <- reconstruct(s, groups = list(trend = c(1, 2)))
  per <- series - eps
  mse <- sum((rec$trend - per)^2)/length(per)
})
teoplitz <- apply(x, 2, FUN = function(series){
  s <- ssa(series, kind = "toeplitz-ssa")
  rec <- reconstruct(s, groups = list(trend = c(1, 2)))
  per <- series - eps
  mse <- sum((rec$trend - per)^2)/length(per)
})

plot(y = teoplitz, x = alpha, type = "l", ylab = "MSE", xlab = 'alpha', main = "Teoplitz vs Basic SSA")
lines(y = basic,  x = alpha,  col = 'red')

legend("topleft", c("Teoplitz", "Basic"), lty=c(1,1), lwd=c(2.5,2.5),col=c("black","red"))

Последовательный анализ реального ряда

uk <- read.csv("UK.csv", sep = ',', stringsAsFactors = FALSE)
uk <- ts(uk$EXPEND, start = c(1980, 1), frequency = 12)
df.uk <- data.frame(seq.Date(from = as.Date("1980-01-01"), to = as.Date("2016-12-01"), by = "month"), uk)
colnames(df.uk) <- c("DATE", "EXP")
head(uk, 24)

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1980 117 116 166 180 202 290 298 441 388 260 175 105
1981 137 142 176 231 240 316 363 537 487 324 185 133

plot(uk)

Периодограмма ряда

spec.pgram(uk, detrend = TRUE, fast = FALSE, log='no', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

Сглаживание

• Тренд сложной формы, возможно, не описывается рядом конечной размерности
• Если выберем большое окно, компоненты тренда смешаются с периодическими
• Возьмем маленькое окно: хорошо выделим тренд, но смешаются остальные компоненты
• Решение: сначала выделим тренд с маленьким окном, потом применим SSA с большим окном к остатку.

Выбираем небольшую длину окна кратную периоду сезонной компоненты, L = 24.

uk.ssa <- ssa(uk, L = 24)

Вклады

plot(uk.ssa)

Видно, что собственные числа периодических компонент близкие. Скорее всего эти компоненты смешаются (проблема с сильной разделимостью).

Собственные вектора

plot(uk.ssa, type="vectors", idx=1:12)

plot(uk.ssa, type="paired", idx=1:11)

Восстановленные по элементарным компонентам

plot(uk.ssa, type ='series', groups = as.list(1:12))

Первая элементарная компонента соответствует тренду, дальше идут периодики.

Выделение тренда

uk.trend <- reconstruct(uk.ssa, groups = list(1))
plot(uk)
lines(uk.trend$F1, col = 'red')

uk.detrended <- residuals(uk.trend)

Периодограмма остатка

spec.pgram(uk.detrended, detrend = FALSE, fast = FALSE, log='no', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

Ожидаем, что дальше будет достаточно выделить 9 компонент (синус с периодом 2 имеет ранг 1).

Сезонность

Возьмем максимальное \(L \leq N/2\) кратное 12, т.е. 216 (длина ряда 444).

uk.detrended.ssa <- ssa(uk.detrended, L=216)
plot(uk.detrended.ssa)

Видно ступеньки, которые, возможно, соответствуют синусоидам.

plot(uk.detrended.ssa, type = "paired", idx = 1:30, plot.contrib = FALSE)

Заметны компоненты, соответствующие периодам 12, 6, 4 и 2.4.

Компонента, соответствующая периоду 2.

plot(uk.detrended.ssa, type = "vectors", idx = 9, plot.contrib = FALSE)

Проверим

parestimate(uk.detrended.ssa, groups = list(1:9), method = "esprit-ls")

   period     rate   |    Mod     Arg  |     Re        Im
    3.995   0.005162 |  1.00518   1.57 | -0.00196   1.00517
   -3.995   0.005162 |  1.00518  -1.57 | -0.00196  -1.00517
   11.991   0.004041 |  1.00405   0.52 |  0.86934   0.50236
  -11.991   0.004041 |  1.00405  -0.52 |  0.86934  -0.50236
    2.000   0.003889 |  1.00390   3.14 | -1.00390   0.00000
    5.995   0.003856 |  1.00386   1.05 |  0.50114   0.86983
   -5.995   0.003856 |  1.00386  -1.05 |  0.50114  -0.86983
    2.401   0.003558 |  1.00356   2.62 | -0.86881   0.50231
   -2.401   0.003558 |  1.00356  -2.62 | -0.86881  -0.50231

W-корреляционная матрица

plot(wcor(uk.detrended.ssa, groups = as.list(1:20)))

uk.seasonality <- reconstruct(uk.detrended.ssa, groups = list(1:9))
res1 <- residuals(uk.seasonality)

par(mfrow=c(4,1))
plot(uk, type = 'l', ylab='Original')
plot(uk.trend$F1, type='l', ylab='Trend')

plot(uk.seasonality$F1, type='l', ylab = 'Seasonal Component')
plot(res1, type='l', ylab='Residuals')

В остатке остался тренд.

acf(res1)

spec.pgram(res1, detrend = FALSE, fast = FALSE, log='no', xaxt = 'n')
axis(1, at = c(0,1,2,3,4,5,6), labels = c('0', '1/12', '2/12', '3/12', '4/12', '5/12', '6/12'))

Оценка дисперсии шума

sigma.n.sqr <- ssa(res1^2, L = 30)
sigma.n <- sqrt(reconstruct(sigma.n.sqr, groups = list(1))$F1)
plot(res1, type='l')
lines(sigma.n, type = 'l', col='blue')

lines(-sigma.n, type='l', col='blue')

Forecasting

Рассмотрим тот же самый ряд, но без последних шести лет.

uk.train <- head(uk, -72)
plot(uk.train)

uk.train.ssa <- ssa(uk.train, L = 24)
uk.train.trend <- reconstruct(uk.train.ssa, groups = 1)
uk.train.detrended <- residuals(uk.train.trend)
uk.train.detrended.ssa <- ssa(uk.train.detrended, L = 180)
uk.train.seasonality <- reconstruct(uk.train.detrended.ssa, groups = list(1:9))

R-forecasting

Продолжаем ряд с помощью min-norm LRR.

uk.trend.forecast <- predict(uk.train.ssa, method = "recurrent", len = 72, groups = 1)
uk.seasonality.forecast <- predict(uk.train.detrended.ssa, method = "recurrent", len = 72, groups = list(1:9))
plot(uk, xlim = c(2000, 2020))
lines(uk.trend.forecast + uk.seasonality.forecast, col = 'red')

MSE

sum((tail(uk, 72) - (uk.trend.forecast + uk.seasonality.forecast))^2)/72

[1] 117129.6

V-forecasting

uk.trend.forecast <- predict(uk.train.ssa, method = "vector", len = 72, groups = 1)
uk.seasonality.forecast <- predict(uk.train.detrended.ssa, method = "vector", len = 72, groups = list(1:9))
plot(uk, xlim = c(2000, 2020))
lines(uk.trend.forecast + uk.seasonality.forecast, col = 'red')

MSE

sum((tail(uk, 72) - (uk.trend.forecast + uk.seasonality.forecast))^2)/72

[1] 76774.87

Доверительные интервалы

Bootstrap:

• Выделяем сигнал с помощью SSA
• Оцениваем дисперсию остатка и много раз моделируем гауссовский шум
• Выделяем сигнал из оцененного сигнала + моделированного шума
• Получим много оценок сигнала
• Строим продолжения
• Отбрасываем \((1 - \gamma)/2\) с концов

uk.trend.forecast <- predict(uk.train.ssa, method = "bootstrap-recurrent", len = 72, groups = 1)
uk.seasonality.forecast <- predict(uk.train.detrended.ssa, method = "bootstrap-recurrent", len = 72, groups = list(1:9))
plot(uk, xlim = c(2010, 2018))
lines(uk.trend.forecast[,1] + uk.seasonality.forecast[,1], col = 'red')

lines(uk.trend.forecast[,2] + uk.seasonality.forecast[,2], col = 'blue',lty="dashed")
lines(uk.trend.forecast[,3] + uk.seasonality.forecast[,3], col = 'blue',lty="dashed")

Заполнение пропусков

Рассмотрим ряд

uk.head <- head(uk, - 120)
plot(uk.head)

Добавим пропуски в ряд

uk.na <- uk.head
missed <- uk.head[121:180]
uk.na[121:180] = NA
plot(uk.na, type = 'l')

Алгоритм

Subspase-based method:

• На этапе декомпозиции, исключаем вектора вложения, содержащие пропущенные элементы
• Из оставшихся векторов составляем траекторную матрицу
• SVD
• Проекция на \(\mathcal{L}_r\)
• Проекция неполных векторов на \(\mathcal{L}_r\).
• “предсказание во внутрь” с помощью лрф, построенной по имеющимся данным. Предсказываем слева и справа с линейно убывающими весами, а потом усредняем
• Диагональное усреднение

Итеративный метод:

• Заполняем пропуски, например, средним
• r и L фиксированы (например выбрали с помощью cv)
• Применим SSA, вставим интересующую нас часть в исходный ряд, снова применим SSA и так далее

Минус: для прогноза не годится.

uk.na.ssa <- ssa(uk.na, L = 72)

Предсказание во внутрь

g <- gapfill(uk.na.ssa, groups = list(1:8))
plot(unclass(uk.head), type = 'l')
lines(y = unclass(g[121:180]), x = 121:180, col = 'red', lty = 2)

MSE

sum((g[121:180] - missed)^2)/60

[1] 10470.2

Итеративный метод

g.it <- igapfill(uk.na.ssa, groups = list(1:8))
plot(unclass(uk.head), type = 'l')
lines(y = unclass(g.it[121:180]), x = 121:180, col = 'red', lty = 2)

MSE

sum((g.it[121:180] - missed)^2)/60

[1] 8747.204

Итеративный метод оказался точнее.

2D-SSA для разложения изображения

Рассмотрим фотографию спутника Сатурна, Дионы (Источник).

Прочитали картинку

dione <- readJPEG("Dione.jpg")
image(dione, col = grey(seq(0, 1, length = 256)))

Добавим шум

set.seed(100)
noise.matrix <- apply(dione, 2, function(x){rnorm(length(x), sd = 0.3)})
dione.noised <- dione + noise.matrix
image(dione.noised, col = grey(seq(0, 1, length = 256)))

2d-ssa

dione.ssa <- ssa(dione.noised, kind = "2d-ssa", L = c(25,25))

Собственные вектора

plot(dione.ssa, type = "vectors", idx = 1:30, cuts = 255, layout = c(10, 2), plot.contrib = FALSE)

w-корреляционная матрица

plot(wcor(dione.ssa, groups = 1:30),scales = list(at = seq(10, 30, 5)))

Сгруппируем компоненты с частотами похожего вида

dione.reconstructed <- reconstruct(dione.ssa, groups = list(I1 = c(1:13), I2 = c(14:15), I3 = c(16:17)))

Накопленные суммы компонент

plot(dione.reconstructed, cuts = 255, layout = c(5, 1), par.strip.text = list(cex = 0.75), type = "cumsum", at = "free")

Shaped SSA

Уберем черный фон

dione.without.black <- apply(dione, 2, FUN = function(col){
  sapply(col, FUN = function(element){
    if(element <= 0.1){
      return(NA)
    } else{
      return(element)
    }
  })
})

image(dione.without.black, col = grey(seq(0, 1, length = 256)))

Добавим шум

dione.without.black.noised <- dione.without.black + noise.matrix
image(dione.without.black.noised, col = grey(seq(0, 1, length = 256)))

Попробуем изменить форму окна на круг.

dione.ssa.shaped <- ssa(dione.without.black.noised, kind = "2d-ssa", wmask = circle(20))

Some field elements were not covered by shaped window. 1253 elements will be ommited

plot(dione.ssa.shaped, type = "vectors", idx = 1:30, cuts = 255, layout = c(10, 2), plot.contrib = FALSE)

w-корреляционная матрица

plot(wcor(dione.ssa.shaped, groups = 1:30),scales = list(at = seq(10, 30, 5)))

По матрице можно сделать вывод о том, что приблизительно 10 первых компонент относятся к сигналу, а остальное это шум.

Сгруппируем компоненты с частотами похожего вида

dione.reconstructed.shaped <- reconstruct(dione.ssa.shaped, groups = list(I1 = c(1:6), I2 = c(7:12), I3 = c(13:19)))

Накопленные суммы компонент

plot(dione.reconstructed.shaped, cuts = 255, layout = c(5, 1), par.strip.text = list(cex = 0.75), type = "cumsum", at = "free")

На вид, исходная картинка лучше отделилась от шума в случае shaped SSA.

Exponential Smoothing

(ES vs SSA vs ARIMA)

es.fit <- es(uk.train, h = 72)

Forming the pool of models based on... ANN, ANA, ANM, AAM, Estimation progress:    45%55%64%73%82%91%100%... Done! 

MSE

sum((tail(uk,72) - es.fit$forecast)^2)/72

[1] 91021.91

Модель

es.fit$formula

[1] "y[t] = (l[t-1] + b[t-1]) * s[t-12] * e[t]"

uk.trend.forecast <- predict(uk.train.ssa, method = "recurrent", len = 72, groups = 1)
uk.seasonality.forecast <- predict(uk.train.detrended.ssa, method = "recurrent", len = 72, groups = list(1:9))
plot(uk, xlim = c(2010, 2018))
lines(uk.trend.forecast + uk.seasonality.forecast, col = 'red')

sum((tail(uk,72) - (uk.trend.forecast + uk.seasonality.forecast))^2)/72

[1] 117129.6

auto.fit.uk <- auto.arima(uk.train, trace = TRUE, stepwise = FALSE, allowdrift = FALSE)

 ARIMA(0,0,0)(0,1,0)[12]                    : 4609.677
 ARIMA(0,0,0)(0,1,1)[12]                    : 4605.92
 ARIMA(0,0,0)(0,1,2)[12]                    : 4570.102
 ARIMA(0,0,0)(1,1,0)[12]                    : 4613.276
 ARIMA(0,0,0)(1,1,1)[12]                    : Inf
 ARIMA(0,0,0)(1,1,2)[12]                    : Inf
 ARIMA(0,0,0)(2,1,0)[12]                    : 4575.118
 ARIMA(0,0,0)(2,1,1)[12]                    : Inf
 ARIMA(0,0,0)(2,1,2)[12]                    : Inf
 ARIMA(0,0,1)(0,1,0)[12]                    : 4507.381
 ARIMA(0,0,1)(0,1,1)[12]                    : 4508.452
 ARIMA(0,0,1)(0,1,2)[12]                    : 4485.141
 ARIMA(0,0,1)(1,1,0)[12]                    : 4519.165
 ARIMA(0,0,1)(1,1,1)[12]                    : 4513.546
 ARIMA(0,0,1)(1,1,2)[12]                    : Inf
 ARIMA(0,0,1)(2,1,0)[12]                    : 4500.653
 ARIMA(0,0,1)(2,1,1)[12]                    : Inf
 ARIMA(0,0,1)(2,1,2)[12]                    : Inf
 ARIMA(0,0,2)(0,1,0)[12]                    : 4437.354
 ARIMA(0,0,2)(0,1,1)[12]                    : 4430.973
 ARIMA(0,0,2)(0,1,2)[12]                    : 4420.724
 ARIMA(0,0,2)(1,1,0)[12]                    : 4438.837
 ARIMA(0,0,2)(1,1,1)[12]                    : 4436.083
 ARIMA(0,0,2)(1,1,2)[12]                    : Inf
 ARIMA(0,0,2)(2,1,0)[12]                    : 4442.367
 ARIMA(0,0,2)(2,1,1)[12]                    : Inf
 ARIMA(0,0,3)(0,1,0)[12]                    : 4420.287
 ARIMA(0,0,3)(0,1,1)[12]                    : 4405.62
 ARIMA(0,0,3)(0,1,2)[12]                    : 4397.304
 ARIMA(0,0,3)(1,1,0)[12]                    : 4411.34
 ARIMA(0,0,3)(1,1,1)[12]                    : 4411.748
 ARIMA(0,0,3)(2,1,0)[12]                    : 4422.04
 ARIMA(0,0,4)(0,1,0)[12]                    : 4401.569
 ARIMA(0,0,4)(0,1,1)[12]                    : 4382.084
 ARIMA(0,0,4)(1,1,0)[12]                    : 4386.684
 ARIMA(0,0,5)(0,1,0)[12]                    : 4383.458
 ARIMA(1,0,0)(0,1,0)[12]                    : 4407.718
 ARIMA(1,0,0)(0,1,1)[12]                    : 4344.179
 ARIMA(1,0,0)(0,1,2)[12]                    : 4340.676
 ARIMA(1,0,0)(1,1,0)[12]                    : 4357.705
 ARIMA(1,0,0)(1,1,1)[12]                    : 4352.702
 ARIMA(1,0,0)(1,1,2)[12]                    : Inf
 ARIMA(1,0,0)(2,1,0)[12]                    : 4366.902
 ARIMA(1,0,0)(2,1,1)[12]                    : 4366.76
 ARIMA(1,0,0)(2,1,2)[12]                    : 4368.57
 ARIMA(1,0,1)(0,1,0)[12]                    : 4345.68
 ARIMA(1,0,1)(0,1,1)[12]                    : 4280.329
 ARIMA(1,0,1)(0,1,2)[12]                    : 4277.132
 ARIMA(1,0,1)(1,1,0)[12]                    : 4290.246
 ARIMA(1,0,1)(1,1,1)[12]                    : 4288.556
 ARIMA(1,0,1)(1,1,2)[12]                    : Inf
 ARIMA(1,0,1)(2,1,0)[12]                    : 4301.595
 ARIMA(1,0,1)(2,1,1)[12]                    : 4302.477
 ARIMA(1,0,2)(0,1,0)[12]                    : 4346.712
 ARIMA(1,0,2)(0,1,1)[12]                    : 4281.757
 ARIMA(1,0,2)(0,1,2)[12]                    : 4278.092
 ARIMA(1,0,2)(1,1,0)[12]                    : 4291.57
 ARIMA(1,0,2)(1,1,1)[12]                    : 4289.544
 ARIMA(1,0,2)(2,1,0)[12]                    : 4302.708
 ARIMA(1,0,3)(0,1,0)[12]                    : 4348.464
 ARIMA(1,0,3)(0,1,1)[12]                    : 4277.812
 ARIMA(1,0,3)(1,1,0)[12]                    : 4291.172
 ARIMA(1,0,4)(0,1,0)[12]                    : 4346.347
 ARIMA(2,0,0)(0,1,0)[12]                    : 4356.483
 ARIMA(2,0,0)(0,1,1)[12]                    : 4294.508
 ARIMA(2,0,0)(0,1,2)[12]                    : 4293.482
 ARIMA(2,0,0)(1,1,0)[12]                    : 4309.562
 ARIMA(2,0,0)(1,1,1)[12]                    : 4304.906
 ARIMA(2,0,0)(1,1,2)[12]                    : Inf
 ARIMA(2,0,0)(2,1,0)[12]                    : 4318.985
 ARIMA(2,0,0)(2,1,1)[12]                    : 4318.889
 ARIMA(2,0,1)(0,1,0)[12]                    : 4347.574
 ARIMA(2,0,1)(0,1,1)[12]                    : 4282.203
 ARIMA(2,0,1)(0,1,2)[12]                    : 4278.389
 ARIMA(2,0,1)(1,1,0)[12]                    : 4292.378
 ARIMA(2,0,1)(1,1,1)[12]                    : 4289.94
 ARIMA(2,0,1)(2,1,0)[12]                    : 4303.347
 ARIMA(2,0,2)(0,1,0)[12]                    : 4347.527
 ARIMA(2,0,2)(0,1,1)[12]                    : 4284.437
 ARIMA(2,0,2)(1,1,0)[12]                    : 4294.897
 ARIMA(2,0,3)(0,1,0)[12]                    : 4349.054
 ARIMA(3,0,0)(0,1,0)[12]                    : 4354.896
 ARIMA(3,0,0)(0,1,1)[12]                    : 4293.657
 ARIMA(3,0,0)(0,1,2)[12]                    : 4291.62
 ARIMA(3,0,0)(1,1,0)[12]                    : 4306.162
 ARIMA(3,0,0)(1,1,1)[12]                    : 4303.185
 ARIMA(3,0,0)(2,1,0)[12]                    : 4316.591
 ARIMA(3,0,1)(0,1,0)[12]                    : 4350.227
 ARIMA(3,0,1)(0,1,1)[12]                    : 4280.178
 ARIMA(3,0,1)(1,1,0)[12]                    : 4293.699
 ARIMA(3,0,2)(0,1,0)[12]                    : 4350.246
 ARIMA(4,0,0)(0,1,0)[12]                    : 4347.801
 ARIMA(4,0,0)(0,1,1)[12]                    : 4286.691
 ARIMA(4,0,0)(1,1,0)[12]                    : 4298.339
 ARIMA(4,0,1)(0,1,0)[12]                    : 4349.489
 ARIMA(5,0,0)(0,1,0)[12]                    : 4350.566

arima.prediction.uk <- predict(auto.fit.uk, 72)
plot(uk, xlim = c(2010, 2018))
lines(arima.prediction.uk$pred, col = 'red')

lines(arima.prediction.uk$pred + arima.prediction.uk$se, lty = 'dashed', col = 'blue')
lines(arima.prediction.uk$pred - arima.prediction.uk$se, lty = 'dashed', col = 'blue')

MSE

sum((tail(uk,72) - arima.prediction.uk$pred)^2)/72

[1] 477969.6